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ABSTRACT 

The pulsar wind nebula associated with PSR J1826— 1334, HESS J1825 — 137, is a bright very high energy 
source with an angular extent of ~ 1° and spatially-resolved spectroscopic TeV measurements. The gamma- 
ray spectral index is observed to soften with increasing distance from the pulsar, likely the result of cooling 
losses as electrons traverse the nebula. We describe analysis of X-ray data of the extended nebula, as well 
as 3-D time-dependent spectral energy distribution modeling, with emphasis on the spatial variations within 
HESS J1825— 137. The multi-wavelength data places significant constraints on electron injection, transport, 
and cooling within the nebula. The large size and high nebular energy budget imply a relatively rapid initial 
pulsar spin period of 13 ± 7 ms and an age of 40 ± 9 kyr. The relative fluxes of each VHE zone can be 
explained by advective particle transport with a radially decreasing velocity profile with v(r) oc r~ - 5 . The 
evolution of the cooling break requires an evolving magnetic field which also decreases radially from the 
pulsar, B(r, t) oc r~ - 7 E{t) x /' 2 . Detection of 10 TeV flux ~ 80 pc from the pulsar requires rapid diffusion of 
high energy particles with r esc m 90 (i?/10 pc) 2 (E c /100 TeV) -1 year, contrary to the common assumption 
of toroidal magnetic fields with strong magnetic confinement. The model predicts a rather uniform Fermi LAT 
surface brightness out to ~ 1° from the pulsar, in good agreement with the recently discovered LAT source 
centered 0.5° southwest of PSR J1826— 1334 with extension 0.6 ± 0.1°. 
Subject headings: pulsars: individual (PSR J1826— 1334) - X-rays: general 



1. INTRODUCTION 

Pulsar wind nebulae (PWNe) confine magnetized pulsar 
winds via a termination shock, whereupon particles are re- 
accelerated to extremely high energies and emit radiation 
across the electromagnetic spectrum. Many PWNe are spa- 
tially resolved in the radio, X-ray, and even very high energy 
(VHE) wavebands, and thereby provide an excellent labora- 
tory to study not only pulsar winds and dynamics, but also 
shock processes, the ambient medium, magnetic field evolu- 
tion, and particle transport. 

An archetypal PWN system is HESS J 1825- 137 associated 
with the energetic pulsar PSR J 1826— 1334, which possesses 
an energy-dependent TeV morphology thought to be caused 
by relic electrons from earlier epochs (Aharonian et al. 2006). 
This system is also a prime example of offset VHE PWNe, 
with the H.E.S.S. flux centered 0.3° south of the pulsar and 
VHE emission extending 1.2° to the south. The TeV emission 
softens from a photon index of T = 1.8 near the pulsar to 
T = 2.5 in the outer reaches of the nebula. 

PSR J1826-1334 was detected in the Jodrell Bank 20 cm 
radio survey (Clifton et al. 1992), with period 101 ms, spin- 
down power E = 2.8 x 10 36 ergs -1 , and characteristic age 
t c = 21.4 kyr. The dispersion measure distance of the pul- 
sar is ~ 3.9kpc (Cordes & Lazio 2002). X-ray observations 
of the PWN with XMM-Newton showed a compact core with 
a hard photon index (r = 1.6io' 2 ) of size 30" embedded in 
a larger diffuse structure of extension ~ 5' extending to the 
south of the pulsar with a softer photon index of T w 2.3 
(Gaensler et al. 2003). Chandra observations by Pavlov et al. 
(2008) revealed the compact core to be composed of a bright 
inner component rj 7" x 3" surrounded by an outer compo- 
nent elongated in the east-west direction. More recent X-ray 
observations with Suzaku revealed diffuse X-ray emission of 
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T ss 2.0 extending 15' south of the pulsar with little sign of 
spectral softening (Uchiyama et al. 2009). HESS J 1825 -137 
has not been detected in the radio (Osborne et al. 2009). 

Using archival VLA observations of PSR J1826-1334 
Pavlov et al. (2008) measured a proper motion of 
440 d,4 km s -1 (with distance d = Ad^ kpc) at a position an- 
gle 100 deg east of north. This places the pulsar birthsite 
some 8' = 9 d^ pc to the west for an age of 20 kyr (see Fig- 
ure 1), and rules out pulsar proper motion as the cause of the 
large southerly extension of HESS J1825— 127. Gaensler et 
al. (2003) advocated that the offset is instead likely the result 
of an asymmetric reverse shock returning from the north, as 
proposed by Blondin et al. (2001) to explain the similar struc- 
ture observed in Vela-X. Indeed, Lemiere et al. (2006) found 
a compelling spatial correlation between a molecular cloud 
discovered in CO studies and the VHE source. The distance 
to this cloud is measured as 4 kpc, and with a total mass of 
«4x 10 5 M & naturally provides the ISM inhomogeneity to 
the north required for the reverse shock scenario. The fact that 
emission extends roughly perpendicular to the pulsar motion 
suggests that the reverse shock velocity must be significantly 
greater than the pulsar proper motion. 

The Fermi Large Area Telescope (LAT) recently detected 
HESS J1825— 137 using 20 months of survey data (Grondin 
et al. 2011). The gamma-ray emission detected by the LAT 
boasts a similar size to the H.E.S.S. source, with extension 
0.56° ± 0.07° for an assumed Gaussian model. The 1 — 100 
GeV LAT spectrum of this source is well described by a 
power-law with a spectral index of 1.4 ± 0.2. Single-zone 
spectral energy distribution (SED) modeling of the X-ray and 
gamma-ray data showed that ss 5 x 10 49 erg injected in the 
form of electrons evolved over 25 kyr could reproduce the 
broadband aggregate spectrum, and that the hard LAT spec- 
trum was consistent with both a simple power-law electron 
injection spectrum, and with a relativistic Maxwellian with a 
power-law tail (Grondin et al. 201 1). 
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In this paper we describe SED modeling of the X-ray and 
gamma-ray data, with particular emphasis on the spatial varia- 
tions in spectra. In Section 2 we overview the X-ray and VHE 
data, Section 3 details the SED model, Section 4 describes re- 
sults of SED fitting, and Section 5 discusses the conclusions 
that can be drawn from model fitting. 

2. DATA ANALYSIS 
2.1. Suzaku Data 

Diffuse Suzaku X-ray emission surrounding PSR 
J 1826- 1334 was analyzed by Uchiyama et al. (2009). 
The 18' Suzaku XIS field of view is broad enough to encom- 
pass the bright core of HESS J 1825 -137, and Uchiyama et 
al. (2009) reported emission nearly to the edge of all four 
XIS chips. Concurrently, a dedicated background pointing 
35' northeast of the initial exposure and outside the H.E.S.S. 
emission region was used for spectral backgrounds. More 
recently, a greater portion of the H.E.S.S. excess was covered 
by a grid of three Suzaku pointings (unpublished at the 
time of submission) to the South and West of the initial 
pointing (see Table 1). A further Suzaku exposure centered 
on PSR J1826— 1334 (also unpublished) expanded the data 
coverage of the PWN core. These five images therefore cover 
k, 40' x 40' of the brightest portion of HESS J 1825 -137, and 
are listed in Table 1 . A number of point sources are observed 
in the Suzaku mosaic, five of which (A-E) were noted by 
Uchiyama et al. (2009). These sources are listed in Table 2, 
as well as displayed in Figure 1. Sources E and N have no 
cataloged counterpart. 

We utilize the standard pipeline screened events, and ana- 
lyze the XIS chips with XSelect version 2.4. The four most 
recent pointings occurred after the loss of the XIS2 chip, so 
these pointings include the two remaining front side illumi- 
nated chips (XIS0, XIS3) as well as the back side illumi- 
nated chip (XIS 1). All observations operate the XIS in normal 
clocking full window mode, with the data split between two 
editing modes: 3x3 and 5 x 5. An exposure corrected mo- 
saic of the five Suzaku exposures is shown in Figure 1 . Diffuse 
emission clearly extends <~ 18' south of the pulsar position, 



TABLE 1 

Suzaku Data 



No. 


ObsID 


Date 


Pointing" 


Location 


Time (ks) 





501045010 


2006-10-19 


276.90, -13.29 


BG 


52.1 


1 


501044010 


2006-10-17 


276.51,-13.71 


H.E.S.S. core 


50.3 


2 


503028010 


2008-10-15 


276.50, -14.01 


Sofobs#l 


57.2 


3 


503029010 


2008-10-17 


276.20, -13.72 


W of obs #1 


57.2 


4 


503030010 


2008-10-19 


276.19,-13.99 


SW of obs #1 


57.2 


5 


503086010 


2009-03-19 


276.57,-13.59 


PSR 


52.1 



a RA, Dec (degrees) 



TABLE 2 
X-ray Point Sources 



Source 


Name/Location" 


Source 


Name/Location" 


A 


2XMM J182557. 9-134755 


H 


2XMM J182539. 9-133130 


B 


2XMM J182617. 1-134111 


I 


2RXP J182535. 2-135049 


C 


2XMM J182629. 5-133648 


J 


2RXP J182458. 3-133747 


D 


2XMM J182620.9-134426 


K 


2RXP J182509. 9-134621 


E 


18:26:12 -13:48:33 


L 


2RXP J182601.4-140423 


F 


2RXPJ182640.4-133911 


M 


2RXP J182436. 5-140429 


G 


2RXPJ182642.6-133625 


N 


18:24:14-13:58:15 



a RA, Dec (Sexagesimal) 



with very little emission beyond this, or to the north of the 
pulsar. 

2.2. H.E.S.S. Data 

The twelve H.E.S.S. spectral extraction regions extend in 
6' chunks 1.2° from the pulsar (0.1° = 7.0 c?4 pc), covering a 
90° sector from approximately 190° — 280°, measured coun- 
terclockwise from North (Aharonian et al. (2006), Figure 4). 
While these regions capture the majority of the VHE emis- 
sion, inspection of the H.E.S.S. excess map indicates that a 
sizable fraction of the emission lies to the southeast of the pul- 
sar and outside of the extraction regions. An accurate value 
of the ratio of X-ray to gamma-ray flux is important for mod- 
eling purposes (see Section 3.1). Therefore, to account for 
the missing VHE flux, and to allow better comparison with 
the Suzaku X-ray data (which extends beyond the 90° sector) 
and LAT GeV gamma-ray data (which is drawn from a very 
large extraction region), we extend the quarter annuli an ex- 
tra 45° such that they now cover 145° — 280°. Barring a full 
reanalysis of the H.E.S.S. data (which is not possible given 
the propriety nature of H.E.S.S. data) the best we can do is 
to scale the amplitude of each slice's spectrum by the ratio 
of total image counts in the original versus extended region. 
As a check of this method, we compute the image counts ra- 
tios between slices. For the original (unsealed) H.E.S.S. re- 
gions, the ratio of image counts between any two regions lies 
within the error bounds of the published H.E.S.S. flux ratios 
between those two regions. This lends some credence to the 
assumption that the flux within a region scales directly with 
the number of counts in the H.E.S.S. excess map. 

The scaling factors are listed in Table 3 (along with X-ray 
scaling factors, described in Section 2.3). These expanded 
regions capture 81% of the counts within the 0.8° radius re- 
gion used to extract the aggregate H.E.S.S. spectrum. The 
VHE flux data points of the SED plots (Figures 6, 7, 9, 11) 
are therefore taken to be the rescaled points from Figure 4 of 
Aharonian et al. (2006), and summing the rescaled flux data 
points of these regions yields data points w 25% below the 
aggregate H.E.S.S. data. This is quite close to the expected 
19% difference expected from the counts ratio in the H.E.S.S. 
excess map. Some of the missing flux is the northeast exten- 
sion of the nebula, but some is likely contributed by adjacent 
regions to the north. Aharonian et al. (2006) estimate the sys- 
tematic error on the flux scale to be 20%, though for SED 
fitting we utilize statistical errors only. Figure 1 shows the 
H.E.S.S. excess map and spectral extraction regions. 



TABLE 3 
H.E.S.S. Spectral Data 



Re 


gion 


X-Ray Scaling H.E.S.S. Scaling 


pa 


Scaled Flux b 


0' 


-6' 


1.13 


1.0 


1.83 ±0.09 


2.9 ±0.3 


6' - 


- 12' 


1.16 


1.5 


1.96 ±0.08 


4.4 ±0.4 


12' 


-18' 


1.27 


1.4 


2.20 ±0.06 


6.0 ±0.4 


18' 


- 24' 




1.4 


2.25 ±0.06 


7.4 ±0.4 


24' 


- 30' 




1.4 


2.32 ± 0.07 


8.7 ± 0.5 


30' 


- 36' 




1.4 


2.37 ±0.06 


9.7 ± 0.5 


36' 


-42' 




1.3 


2.36 ±0.08 


7.4 ±0.5 


42' 


-48' 




1.4 


2.41 ±0.09 


7.6 ±0.6 


48' 


-54' 




1.2 


2.42 ± 0.09 


6.1 ± 0.5 


54' 


-60' 




1.2 


2.59 ±0.13 


6.2 ±0.6 


60' 


-66' 




1.0 


2.43 ±0.09 


5.9 ±0.5 


66' 


- 72' 




1.0 


2.45 ± 0.35 


2.9 ±0.7 



a Data taken from Aharonian et al. (2006), Table 2 

b 0.25 - lOTeV fluxes in units of 10 -12 erg cm" 2 s" 1 
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FIG. 1 . — Left: Suzaku exposure corrected 1 — 8 keV counts mosaic of the front side illuminated chips for all five pointings, smoothed with a 1.2' Gaussian. 
Point sources (see Table 2) are denoted by 1.8' radius dashed circles, and the three X-ray spectral extraction regions are shown in solid white. PSR J1826— 1334 
is indicated by the black cross, and the green wedge indicates the range of possible pulsar birthsites (Pavlov et al. 2008). This green wedge extends some 16', 
corresponding to the presumed birthsite for an age of 40 kyr. The location of each of the five exposures is indicated by the map in the upper right. Right: H.E.S.S. 
excess counts map showing the VHE spectral extraction regions (solid white), and larger extrapolated regions (dashed white). The Suzaku field of view is shown 
in cyan. The bright point source to the south is the 7-ray binary LS 5039. 



2.3. X-ray Spectrum 

Motivated by SED modeling, we extract Suzaku X-ray 
spectra from 6' wedges to mimic the regions of Aharonian 
et al. (2006). Spectral fitting is accomplished with Xspec ver- 
sion 12.6. For background subtraction we follow the same 
procedure as Uchiyama et al. (2009) (a more in-depth discus- 
sion can be found in this paper) and define the same region as 
these authors: a 9' circle on the background exposure, excis- 
ing 1.8' radius circles around the four point sources noted by 
these authors. The non-uniformity of the background is taken 
into account internally within Xspec. Spectra are extracted 
with XSelect, while ARE and RMF response files are gen- 
erated with the xisrmfgen and xissimarfgen functions, which 
internally take into account the different responses of each 
chip. The ASCA FTOOL addascaspec is utilized to com- 
bine the spectra and responses of the XIS front side illumi- 
nated chips. The XIS1 back side illuminated chip possesses a 
markedly different response function from the front side illu- 
minated chips, and so is analyzed in parallel rather than added 
to the other chips. 

Initially, we extract and fit the spectrum of the 1.5' circle 
surrounding PSR J 1826— 1334 titled Region A by Uchiyama 
et al. (2009). Fitting an absorbed power-law model we 
achieve a fit well within errors of the parameters reported by 
these authors (r — 1.7 ± 0.1, unabsorbed 1 — 8 keV flux 
of 1.2 ± 0.2 x 10~ 12 erg cm -2 s _1 ). We subsequently de- 
fine regions mimicking the inner 6' circle surrounded by 6' 
thick partial annuli shown in Figure 1. From these regions 
we excise 1.8' radius circles around the 5 point sources de- 
fined by Uchiyama et al. (2009). For the innermost 6' circle 
we use both the original dataset (501044010), as well as the 
most recent dataset (503086010) and jointly fit the absorption 
and power-law index. Since the original dataset does not en- 
compass the entirety of the 6' circle we allow the power-law 



normalization to vary between the two datasets, though we 
quote the flux from the entire circle. The 6' — 12' region fol- 
lows the same procedure, using both datasets 501044010 and 
503086010. For overlapping regions (exposure 1 and expo- 
sure 5) the joint spectral fits to both exposures are completely 
consistent (albeit with smaller errors) in both index and to- 
tal flux (when taking into account the differing areas of the 
regions) with fits to the individual exposures. 

All fits are to an absorbed power-law from 1 — 8 keV with 
quoted single parameter 68% (ler) errors. We fit the 6' — 12' 
region with Nh fixed at the value determined in the inner- 
most region, and find values completely consistent (albeit 
with smaller errors) with fitting with Nh free. Though Nh 
may vary over the nebula, we see no evidence of this, and 
IR maps show no appreciable gradient over the nebula, so for 
simplicity we adopt a uniform column density. The pulsar will 
contribute flux to the innermost region, though with a 1 — 8 
keV flux of - 2 x 10~ 14 erg cm- 2 s" 1 (Pavlov et al. 2008), 
the pulsar presents only a minor perturbation to the PWN flux. 

The 12' — 18' wedge stretches across three different point- 
ings and suffers from a much lower surface brightness than the 
inner two regions, greatly complicating flux extraction from 
this region. Due to the limited statistics, we fix the absorption 
to the best fit value from the innermost region. Extracting flux 
from the portion of this region which lies within the original 
pointing (501044010) yields a ~ 5<r detection. No signifi- 
cant emission further south and west of the pulsar is seen in 
any of the three exposures. In fact, we can use the apparent 
lack of flux in the southwest exposure furthest from the pulsar 
(503030010) to estimate the uncertainty in the background. 
Selecting a background region from this chip alters flux esti- 
mates by « 1% for Region 1 (0' - 6'), « 3% for Region 2, 
and pa 26% for Region 3. An additional source of error stems 
from the uncertainty in the ARF calculation for complex ex- 
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TABLE 4 
X-Ray Spectral Fits 



Region 




r 


Unabs. Flux b 


X 2 /d.o.f. 


0' -6' 


1.26 ± 0.06 


1.93 ±0.05 


4 an+ 034 
^■ yu -0.32 


196/249 


6' - 12' 
6' - 12' 


1 12+ - 11 
1.3 C 


2.03 ±0.09 
2.17 ±0.05 


3 0'5+ - 39 
3.18 ±0.15 


144/173 
146/174 


12' - 18' 


1.3 C 


1.78 ±0.18 


1 04+ ' 22 


44/59 



a interstellar absorption X 10 22 cm 2 

b 1 — 8keV fluxes in units of 10~ 12 crgcm~ 2 s" 1 

c held fixed 
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FIG. 2. — Radial flux profile with la statistical error bars. Suzaku 1 — 8 
keV X-ray flux is shown in red, while the far greater extent of the H.E.S.S. 
0.25 - 10 TeV flux (from Table 3) is in blue. 

tended sources, which leads to uncertainty in the source flux 
(Ishisaki et al. 2007), which we estimate at m 10%. Note that 
background subtraction errors dominate for Zone 3; in fact, 
the relatively hard spectrum fit here may reflect residual back- 
ground contamination. In the SED plots we show regrouped 
X-ray data with the ARF and background systematic errors 
summed in quadrature with la statistical errors. 

To be perfectly consistent, for modeling purposes we 
should utilize identical regions for both X-ray and VHE 
gamma-ray data. The excised point sources and slight por- 
tion of the H.E.S.S. spectral extraction regions which extends 
beyond the Suzaku exposures therefore require a small adjust- 
ment to the X-ray flux level. The raw X-ray data is therefore 
rescaled by the fraction of missing area, with the scaling fac- 
tors shown in Table 3. Table 4 and the SED plots therefore 
display the rescaled fluxes (with errors adjusted accordingly) 
in order to remain fully consistent with the H.E.S.S. data. 

The X-ray flux falls off rapidly with radial distance from the 
pulsar, likely the result of both lower magnetic fields as well 
as cooling losses, as discussed below. Figure 2 shows a radial 
profile of fluxes in the X-ray and TeV energy bands, clearly in- 
dicating the much greater extent of the nebula in the H.E.S.S. 
regime. For this figure, we fix Nh at 1.3 x 10 22 cm -2 and 
extract X-ray spectra from 3' width regions for greater detail 
than the 6' regions shown in Table 4. 



3. SED MODELING 

3.1. Motivation 

Grondin et al. (201 1) applied a one-zone SED model to the 
new Fermi LAT data of HESS J1825— 137 and found an ad- 
equate fit to the aggregate gamma-ray and X-ray data. One- 
zone models can help constrain certain global nebula prop- 
erties, such as the total energy injected in the form of elec- 
trons, the electron spectrum in the uncooled regime, and the 
mean magnetic field. Yet such models by definition assume a 
completely homogeneous nebula, contrary to the differing X- 
ray and VHE sizes and spectral steepening observed in both 
X-rays and gamma-rays. To further investigate these phe- 
nomena we implement a spatially resolved model of HESS 
J1825-137. 

We begin by investigating how simple scaling relations in- 
form model implementation. Starting assumptions include: 
the source of all particles and magnetic field in the nebula is 
PSR J1826— 1334, ambient photon fields are static, and the 
nebula is 4 kpc distant. Let us first consider the electrons 
present in the innermost zone (Zone 1) adjacent the pulsar. In 
a transverse magnetic field these particles synchrotron radiate 
photons at mean energy: 

E 1 « 2.2(£ e /100TeV) 2 G B /10^G) keV - (!) 

We scale with a low field of 10 fiG since, as we shall see later, 
even in the innermost zone the mean field is only ~ 10 /iG for 
the best fit model. Such particles cool rapidly, with a lifetime 

(t = E/E) of: 

r sync a 820 (tfe/lOOTeV^OB/lO/xG)- 2 year. (2) 

Therefore the cooling timescale for electrons radiating syn- 
chrotron photons at mean energy is; 

r sync k 1200( J B 7 /lkeV)" 1/2 (S/10^G)" 3/2 year. (3) 

The flat (T ps 2) X-ray spectrum extending to 10 keV of 
Region 1 indicates that the electrons responsible for this 
emission are largely uncooled, implying r < T sync ~ 
380 (B/10 /iG)~ 3 ' 2 year. The inner region therefore must 
be dominated by particles recently injected by the pulsar. In 
the uncooled regime the synchrotron spectral flux scales with 
the magnetic field and number of electrons as: F(E) sync oc 
ri e x B z / 2 . Over timescales short compared to the 21 kyr spin- 
down age of the pulsar the number of electrons n e in Zone 1 is 
roughly proportional to the time period over which electrons 
have been injected St, so F(E) sync oc St x £? 3 / 2 , a relation 
we will return to momentarily. 

Electrons inverse Compton (IC) upscatter CMB photons to 
a mean energy of: 

E 1 « 0.32 (£ e /10 TeV) 2 TeV. (4) 

For young electrons IC cooling is of little consequence, since 
10 TeV electrons scattering off the CMB cool in tic ~ 
140 (E e /10 TeV)' 1 kyr. We scale to 10 TeV rather than 100 
TeV since at 100 TeV the Thomson regime is no longer real- 
ized and Klein-Nishina corrections must be incorporated. For 
example, while the above scaling relation for inverse Comp- 
ton cooling holds for electrons with E e < 10 TeV, 100 TeV 
electrons cool slowly over 36 kyr, rather than the 14 kyr 
predicted by this relation. We therefore consider only syn- 
chrotron cooling in the following discussion. The lowest en- 
ergy H.E.S.S. point at E 1 w 0.3 TeV requires electrons of 
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energy E e w 10 TeV, which synchrotron cool slowly at a rate 
of T sync w 8200 (5/10 /xG) -2 year. Therefore, electrons re- 
sponsible for the low energy H.E.S.S. points cool quite slowly 
compared to those responsible for the X-ray flux. 

The electrons responsible for the lowest energy H.E.S.S. 
points (and therefore all the Fermi LAT flux points as well) 
are radiatively uncooled over time spans of a few thousand 
years, and provide a good indication of just how many elec- 
trons are present in the inner nebula. This is due to the fact 
that for E e < 10 TeV the Thomson limit holds and IC flux 
scales linearly with the number of electrons (n e ) present: 
F(E)jc oc n e . The exact number of electrons at 10 TeV 
will depend on the electron injection spectrum parameters as 
well as the injection time St, and very weakly on the magnetic 
field due to feeble synchrotron cooling, so F(E)jc oc St. 
The X-ray data breaks the degeneracy in magnetic field since 
F(E) sync oc St x £? 3 / 2 , so between X-ray and gamma-ray 
data the electron injection spectrum, age, and magnetic field 
are very well constrained. For modeling the data, this means 
that once an injection spectrum is selected, the gamma-ray 
flux determines the injection time, while the ratio of X-ray to 
gamma-ray flux determines the magnetic field. 

3.2. Pulsar Spin-down 

As pulsars spin down, they dissipate rotational kinetic en- 
ergy according to: 



E B = {l-r])E. 



(10) 



e = mn, 



(5) 



with £1 the angular frequency and / the neutron star's moment 
of inertia, assumed to be 10 45 gem 2 . This energy goes into a 
magnetized particle wind such that for a braking index of n: 



n oc fi™. 



(6) 



For spin-down via magnetic dipole radiation n = 3, though n 
has only been confidently measured for five pulsars, in each 
case falling between 2 < n < 3 (Livingstone et al. (2007) and 
references therein). Integrating Equation 6 yields the age of 
the system (Manchester & Taylor 1977): 



P 



(n-l)\P 



1- ^ 



(7) 



where P is the initial spin period, and P the period deriva- 
tive. For P <C P and n = 3 this equation reduces to the 
characteristic age of the pulsar r c = P/2P. The spin-down 
luminosity of the pulsar evolves as (Pacini & Salvati 1973): 



E(t)=E (l + ^j < ""\ 
with the initial spin-down timescale defined as: 



TO = 



(n - 1) | Po 



(8) 



(9) 



Given that the current P, P, and E are known, once an initial 
period P and braking index n are selected the age T and 
spin-down history of the system are determined according to 
the equations above. 

The energy flux in the magnetized wind E is split between 
electromagnetic Eb and particle energy E e , with: 

E e = r)E 



Their ratio is commonly referred to as the wind magnetization 
parameter: 



a = E B /E e = (1 - r,)/r], 



(11) 



with a typically thought to be <C 1 past the termination 
shock. The particle energy is presumed to go into an elec- 
tron/positron plasma with an exponentially cutoff power-law 
spectrum: 



% oc ET>e-*l*~* 
dE 



E > E min , (12) 



with E min the minimum particle energy. We initially select 
Emin as 100 GeV, which is well within the realm of minimum 
particle energies considered by Kennel & Coroniti (1984b) in 
magnetohydrodynamical modeling of the Crab Nebula. The 
normalization of this power-law varies with time following 
Equation 8, though we treat the index, minimum energy, and 
cutoff energy as static. We begin by modeling pure power-law 
injection, although we later consider an injection spectrum 
composed of a relativistic Maxwellian with a non-thermal tail, 
as proposed by Spitkovsky (2008). 

3.3. Spatial Evolution 

Pulsar winds flow out at relativistic speeds until the ram 
pressure of the wind is balanced by the internal pressure of the 
PWN at the termination shock. High-resolution X-ray images 
can help illuminate the location of this termination shock, and 
Chandra images of the pulsar reveal that the brightest inner 
PWN component stretches some 7" x 3" (0.14 x 0.06 pc 2 at 
4 kpc) (Pavlov et al. 2008), with the short axis roughly ori- 
ented toward the pulsar birthsite. The overall extent of the 
equatorial flow is generally 2 — 3 times larger than the termi- 
nation shock radius (Ng & Romani 2004), and we accordingly 
adopt a value of r ts = 0.03 pc for the radius of the termination 
shock. 

Past the termination shock the pulsar wind expands radially 
at a much reduced velocity. In the model of of Kennel & Coro- 
niti (1984a) the rapid rotation of the pulsar wraps up magnetic 
field lines, resulting in a primarily toroidal field at and beyond 
the termination shock. The field direction is therefore approx- 
imately perpendicular to the wind flow velocity. Since parti- 
cles effectively cannot cross toroidal magnetic field lines, the 
bulk flow of particles predicts an ordered layering of particle 
ages with increasing distance from the pulsar, of age: 



t(r) 



vir'^dr'. 



(13) 



with v(r) the velocity profile. 

On larger scales, we consider the evolution of the parent 
supernova remnant (SNR) of PSR J1826-1334. Early in 
the pulsar lifetime the SNR expands freely into the ambient 
medium, typically at a velocity of <~ 10, 000 km s -1 . The 
large ~ 80 pc size of the VHE PWN and lack of any ob- 
served SNR shell led de Jager & Djannati-Atai (2008) to pro- 
pose a very large SNR of radius w 120 pc resulting from 
an energetic supernova explosion (Esn = 3 x 10 51 erg) of 
age w 40 kyr, and a very low ambient medium density of 
n w 0.001 cm~ 3 . Early on, a spherically symmetric PWN 
expands rapidly into unshocked supernova ejecta with nearly 
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constant velocity (van der Swaluw et al. 2001): 



E 



1/5 



'V» „ --••'.100 . , 

\ I0 iy ergs -1 / \3xl0 51 erg 



Esn 



3/10 



8M, 



1 kyr 



As the SNR expands and sweeps up more mass, it even- 
tually reaches a point where the swept up ISM material is 
approximately equal in mass to the mass of the supernova 
ejecta. At this point the SNR enters the Sedov-Taylor phase 
as the swept-up material begins to dominate the SNR dynam- 
ics. Adopting a uniform density ISM, for the standard core- 
collapse supernovae assumption of a constant density inner 
core surrounded by a p oc r~ 9 outer envelope the transition to 
the Sedov-Taylor phase occurs at (Truelove & McKee 1999): 



ts- 



:7100 



Esn 



-1/2 



3 x 10 51 erg 
n \-V3 



M, 



8M m 



5/6 



vO.OOl cm- 



year. 



(15) 



In this phase a reverse shock is driven deep into the SNR in- 
terior, and for a SNR expanding into a homogenous medium 
the reverse shock will impact the PWN, crushing it and lead- 
ing to a series of reverberations which can last many thou- 
sands of years, as demonstrated compellingly in Gelfand et 
al. (2009). Calculations by van der Swaluw et al. (2004) 
yielded an upper limit on the reverse shock collision with the 
PWN at t rs < 5 ts-T- A supernova explosion into an evacu- 
ated cavity of no = 0.001 cm -3 therefore yields a maximum 
PWN-reverse shock collision time of w 35 kyr (nearly double 
the characteristic age of the pulsar), implying that the PWN 
could to this day be rapidly expanding to the southwest into 
unshocked SNR ejecta. Of course, reflection off of the molec- 
ular cloud discovered by Lemiere et al. (2006) could result in 
a reverse shock collision on a timescale reduced by an order 
of magnitude from this estimate. 

Axisymmetric 2-D simulations by Blondin et al. (2001) 
showed that SNR expansion into an inhomogeneous medium 
leads to a PWN offset from the pulsar position and mixing of 
the thermal gas of the supernova ejecta and the relativistic gas 
of the PWN. Rayleigh-Taylor instabilities aid in the mixing 
and disruption of the nebula as bubbles of swept up material 
penetrate the PWN. The simulations of Blondin et al. (2001) 
showed a cometary nebula ~ 20 kyr after the reverse shock 
interaction, and nearly complete disruption of the nebula after 
<~ 50 kyr. These simulations do not take into account the mag- 
netic field that is mixed in with the relativistic particles, and 
which can dampen instabilities along magnetic field lines. In- 
deed, recent simulations by Bucciantini et al. (2004) suggest 
that the growth of instabilities is highly suppressed for high 
a flows, and significant structure may exist even some <~ 20 
kyr after reverse shock interaction. This suggests that the or- 
dered layering of particle ages predicted by Equation 13 may 
remain valid in the case of an anisotropic reverse shock. 

If one invokes PWN expansion within a freely expand- 
ing SNR to account for the ~ 80 pc size of the PWN, a 
mean velocity of 3100 km s -1 over 25 kyr is required. While 
this velocity is quite high for a PWN, it is permitted by 
Equation 14 for a very low density ambient medium. The 
17' (20^4 pc) displacement of the H.E.S.S. centroid south 



of PSR J1826— 1334 (Aharonian et al. 2006) cannot be ex- 
plained by pure expansion, however, and requires some mech- 
anism to achieve this offset. 

One possibility for the southerly displacement is the exis- 
tence of a powerful jet rapidly transporting particles far from 
the pulsar. X-ray images show no indication of any such fea- 
ture, however. In addition, PWN jets are typically aligned 
with the pulsar spin axis, and a southern jet would be roughly 
perpendicular to the inferred pulsar proper motion, also typi- 
cally aligned with the pulsar spin axis. 

A more palatable alternative is that the southerly offset from 
the pulsar is due to crushing by the reverse shock. In this pic- 
ture the reverse shock quickly reflects off of the molecular 
cloud to the north, and then impacts the PWN in a few thou- 
sand years. The mechanism responsible for the PWN dis- 
placement, be it direct interaction with the reverse shock or a 
resultant anisotropy from the shock, is evidently still at work 
since the young X-ray emitting electrons are clearly concen- 
trated south of the pulsar. For a reverse shock returning from 
the North, particles initially north of the pulsar will be dis- 
placed south and mixed throughout the nebula. The fact that 
the H.E.S.S. spectral index of HESS J1825-137 increases 
monotonically from the pulsar hints that such mixing is likely 
far from complete. Precisely how the reverse shock crush- 
ing and displacement alters any ordered layering of the PWN 
requires detailed 2-D or 3-D magnetohydrodynamical simu- 
lations and is beyond the scope of this paper. We ignore such 
effects, and hence tacitly assume that for HESS J1825— 137 
reverse shock interactions serve primarily to displace rather 
than compress and mix the nebula. 

3.4. Morphology 

The shape of the nebula as viewed from earth resembles 
a circular sector. While the 3 -dimensional morphology of 
HESS J1825 — 137 is uncertain, if the morphology results 
from a shock returning from the north, this implies nebular 
axial symmetry about the normal to the shock front. One sim- 
ple geometry which matches the Suzaku and H.E.S.S. mor- 
phology is a spherical wedge of opening angle 135° (solid 
angle 3.88 steradians) with symmetry axis along the H.E.S.S. 
symmetry axis 17° counterclockwise from North. The molec- 
ular cloud which presumably gives rise to the assymetric re- 
verse shock is located at the same 4 kpc distance as the pul- 
sar, so we assume that the shock normal is perpendicular to 
the line of sight. For modeling purposes we split this wedge 
into twelve 3-D zones which mimic the spectral extraction re- 
gions. Given the axisymmetry, the twelve spatial zones form 
a series of nested bowls, or partial spherical wedges, of thick- 
ness 7.0d,4 pc. Similar to the spectral extraction regions, we 
label Zone 1 adjacent to the pulsar (extending — 7 pc), and 
Zone 12 the outermost spherical wedge (77 — 84 pc from the 
pulsar). Henceforth we refer to Regions as the 2-D projections 
on the sky of the 3-D spatial Zones. Figure 3 demonstrates the 
adopted morphology, with the symmetry axis oriented verti- 
cally for clarity. The final volume of each bubble is deter- 
mined by the spherical wedge morphology discussed above. 
These volumes are listed in Table 5. 

When projected into the sky, the counts from each bowl 
will overlap into multiple spectral extraction regions. In or- 
der to determining the degree of overlap we perform a Monte 
Carlo integration by projecting the 3-D spatial zones onto the 
plane of the sky. For example, consider Zone 8 extending 
49 pc — 56 pc from the pulsar. We populate this zone with a 
large number of random points, then project each point onto 
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FIG. 3.— Top: Spatial model of HESS J1825-137 as seen from Earth, 
with the symmetry axis oriented vertically. The pulsar is located at the origin, 
and the axes indicate the extent of nebula in parsecs. Bottom Left: One of 
the twelve spatial zones (Zone 8: 49 - 56 pc from the pulsar) plotted with 
a slightly elevated viewpoint to demonstrate the concave nature of the zone. 
Bottom Right: The same spatial zone seen edge-on from the Earth line of 
sight. 




FIG. 4. — Smoothed projected counts from Zone 8 extending 49 — 56 pc 
from the pulsar. The 6' width spectral extraction regions are overlaid, demon- 
strating the overlap of counts on the interior regions. The location of the 
pulsar is marked with a cross. 

the appropriate plane defined by the viewing angle. We then 
overlay the spectral extraction regions onto this projection and 
determine the fraction of counts which lie within each bin. 
Assuming each bin is spatially uniform, the fraction of counts 
within each region is equivalent to the fraction of flux. For 
all zones save the innermost, significant flux from each spa- 
tial zone is found in the interior spectral extraction regions, as 
demonstrated by Figure 4. Table 5 shows the fraction of flux 
from each zone which falls within each spectral region. 

3.5. Spatial Model 

Motivated by the morphology adopted in the previous sec- 
tion, we treat the PWN as a series of twelve expanding bub- 
bles with the radial extent of each bubble accounted for by 
a preferential southerly expansion. As per the discussion in 
Section 3.3, this southerly expansion is assumed to result from 
an asymmetric reverse shock returning from the north. Parti- 
cles are injected into each bubble for a period of time 5ti, after 
which the bubble is assumed to detach from the pulsar and di- 
rect injection of particles is discontinued. Given the distinct 
possibility that in the southerly direction SNR is still freely 
expanding (leading to nearly free expansion for the PWN), 
for simplicity we assume constant velocity expansion of the 
outer boundary of the nebula v ou t er — 84pc/T, with T the 
pulsar age determined by Equation 7. Despite the adoption 



of a constant outer expansion velocity, particles will still de- 
celerate over time as they traverse the nebula and encounter 
supernova ejecta. In the interior of the nebula we therefore 
consider a radial velocity profile 



v(r,t) 



(16) 



with a a fit parameter, and the time dependence determined 
by the requirement that v ou ter remain constant. Integrating 
Equation 16 yields the position of various bubbles with time. 
The final extent of each bubble (e.g. Zone 1 covers — 7 
pc from the pulsar) determines the required injection time for 
each bubble via Equation 13. 

Once injection ends, the bulk evolution of each bubble is 
primarily adiabatic, as IC cooling is relatively weak and the 
magnetic field is reduced to the point that synchrotron losses 
are severe only for the highest energy particles. For relativis- 
tic particles with adiabatic index 4/3 the particle pressure is 
Ppwn = E e /(3V pwn ), and for purely adiabatic expansion 

4/3 

PpwnVpwn = constant. For all time steps particle adiabatic 
losses are therefore computed via 



E, 



oc V~ 1/3 . 



(17) 



The sound speed in the relativistic gas of the PWN (c/y/S) 
is far greater than the expansion velocity of the PWN, im- 
plying that pressure differences will rapidly equilibrate and 
that the PWN is nearly isobaric. The energy lost due to adia- 
batic expansion goes into accelerating supernova ejecta at the 
PWN-SNR interface, and so we consider adiabatic losses to 
be shared equally among all particles present in the nebula, 
with the fraction of energy lost by each particle a function 
only of the volume evolution of the aggregate nebula. 

3.6. Magnetic Field Model 

In principle, one can use the injected magnetic field en- 
ergy (Equation 10) and spatial evolution of the PWN zones 

to compute the magnetic field (Eb = ^ V, see Gelfand et al. 
(2009)). Due to the uncertainty in how the magnetic field is 
distributed and transported, and the possibility of field ampli- 
fication in the post-shock flow, we adopt a somewhat simpler 
magnetic field evolution. 

A static and uniform magnetic field presents the simplest 
possible magnetic field structure, yet the nebular magnetic 
field is expected to decrease from the termination shock to 
the far reaches of the nebula, as well as decrease with time 
as the pulsar spins down. Near the pulsar one common sim- 
ple formulation is a dipolar magnetic field (B oc r~ 3 ) out to 
the pulsar light cylinder, with a toroidal field (E> cx r _1 ) out to 
the termination shock. The surface magnetic field is estimated 
from the period and period derivative as 



B s =3.2 x 10 19 (PP) 1/2 G. 



(18) 



For PSR J1826-1334 this gives a current value of B s = 2.8x 
10 12 G for a surface radius r s — 10 km. The light cylinder 
is currently located at r; c = c/Sl = 4.8 x 10 8 cm. At the 
wind termination shock the field is amplified by a factor of 
~ 3, which gives B ts ~ 400 [iG for r ts — 0.03 pc, and the 
assumed geometry of B cx r~ 3 out to the light cylinder and 
B <x r _1 from the light cylinder to the termination shock. 
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TABLE 5 

Projection of Spatial Zones onto Spectral Extraction Regions 





Reg 1 


Reg 2 


Reg 3 


Reg 4 


Reg 5 


Reg 6 


Reg 7 


Reg 8 


Reg 9 


Reg 10 


Reg 11 


Reg 12 V (pc 3 ) 


Zone 1 


1.00 






















440 


Zone 2 


0.14 


0.86 




















3100 


Zone 3 


0.01 


0.30 


0.69 


















4800 


Zone 4 




0.07 


0.34 


0.59 
















16000 


Zone 5 




0.01 


0.13 


0.34 


0.52 














27000 


Zone 6 






0.05 


0.15 


0.32 


0.47 












40000 


Zone 7 






0.02 


0.08 


0.16 


0.30 


0.44 










56000 


Zone 8 








0.04 


0.09 


0.16 


0.30 


0.41 








75000 


Zone 9 








0.02 


0.06 


0.10 


0.15 


0.28 


0.39 






96000 


Zone 10 








0.01 


0.04 


0.07 


0.10 


0.15 


0.27 


0.37 




120000 


Zone 1 1 










0.02 


0.05 


0.07 


0.10 


0.15 


0.26 


0.35 


150000 


Zone 12 










0.01 


0.03 


0.05 


0.07 


0.10 


0.15 


0.25 


0.34 180000 



NOTE. — Each row gives the fraction of flux from each spatial zone falling within each spectral region. For example, Column 1 , Row 2 indicates the fraction of flux from Zone 2 
which falls within Region 1 (0.141). The flux within Region N is the sum of the fluxes in zones 1-12 times the coefficients (C^/v) in column N, Fjv — Y2i=i Ci,N F%- 



The termination shock occurs at radius: 

- (sifc) * ■ <19) 

where e is the filling factor of the pulsar wind. For an isotropic 
wind e = 1, though more common is a scaling with co- 
latitude 9 of sin 2 6, yielding e = 3/2 (Bogovalov & Khang- 
oulyan 2002). The key feature here is that the termination 
shock radius scales as r ts oc (E j P. pwn Y^ 2 . One can esti- 
mate the PWN pressure as P pwn oc p e jV 2 (Equation 26), with 
p e j oc t~ 3 (Equation 24). The assumed constant velocity ex- 
pansion therefore gives P pwn oc t~ 3 , which is identical to 
the Pp Wn oc t~ 30 scaling computed by the Gelfand et al. 
(2009) model prior to reverse shock interactions. The ter- 
mination shock radius accordingly scales roughly as r ts oc 
-E^i) 1 / 2 ^/ 2 , and for a braking index near n = 2 (where 
E(t) oc t~ 3 ) the termination shock remains approximately 
static. Furthermore, since our SED model focuses on emis- 
sion from the entire nebula (which dwarfs r ts ), the precise 
value of r ts has little effect on the resultant SED. For sim- 
plicity, we therefore assume a constant value of r ts = 0.03 
pc. 

Adopting a dipolar followed by toroidal magnetic field ge- 
ometry out to the termination shock, B ts = 3Bi c r\ c r^ 1 = 
3B S r 3 r~ 2 r^ 1 oc p- 3 ' 2 P 1 ' 2 oc E(t)^ 2 . Therefore B ts oc 

■E(i) 1//2 was much higher at earlier times when the pulsar 
spin-down was greatest. Past the termination shock the field 
continues to decrease, though likely at a slower rate than 
B oc 1/r, since this would put the field at < 1 /iG in the 
outermost reaches of the nebula. Instead, we adopt a behavior 
of: 

B(r,t) = 400 ( — r — \ [ - r^rUi I (20) 
v ; V0.03pc/ v 8 x / 

past the termination shock, where f3 is a fit parameter and the 
400 fiG coefficient is the approximate value at the termination 
shock. We should emphasize that this magnetic field model 
does not conserve magnetic flux, since we assume that some 
magnetic field may be generated in the post-shock flow. The 
mean magnetic field within the boundaries of each zone is 
determined via a volume integration of Equation 20, and this 
is the value used at each time step since each zone is treated 



as homogeneous in the SED code. 

3.7. Diffusion 

A certain fraction of the injected particles will be lost from 
the nebula due to diffusive escape. For each zone we compute 
the diffusive flux between zones by solving: 

J i (E e ,t) = -D dn f^ t) cm- 2 s-\ (21) 

with Ji(E e ,t) the diffusive flux in particles per cm 2 s _1 for 
zone i, rn(E e , t), the particle density, and D the diffusion co- 
efficient. The concentration gradient drii (E e ,t) /OR is evalu- 
ated as the concentration difference between adjacent zones 
divided by the distance between zone midpoints. Solving 
Equation 21 for each zone at every time step gives the num- 
ber of particles diffusing between zones. Hence, time permit- 
ting particles typically diffuse sequentially from inner zones 
to outer zones, stepping from Zone 1 to Zone 12. The PWN- 
SNR interface will likely trap the majority of high-energy par- 
ticles within the PWN. In addition, the surrounding SNR is 
expected to contain significant numbers of high energy par- 
ticles. We therefore expect that particles will diffuse rather 
slowly out of the PWN and into the SNR, and so we arbitrar- 
ily set the density in the surrounding SNR nsNR(E e ,t) — 
0.8 x n\i{E e , t), which allows particles to slowly escape the 
PWN from Zone 12. The final SED depends very weakly on 
nsNR.{E ei t), however, given the high degree of cooling and 
large final size of the nebula. 

Toroidal magnetic field lines tend to contain particles very 
effectively, and even with relatively rapid Bohm diffusion 
r esc w 40,000(i?/10pc) 2 (E c /100TeV)- 1 (B/10/xG) year 
(de Jager & Djannati-Atal 2008). In the Bohm limit the cross- 
field mean free path A equals the particle gyroradius, which 
is exceedingly small: A = 0.01(^/100 TeV) (B/IO^G)- 1 
parsec. If the field line structures contain radial components, 
however, particles diffuse far more rapidly. Turbulence and 
mixing caused by the passage of the reverse shock might pro- 
vide the necessary disruption to the magnetic field structure to 
achieve radial field components. Even in young PWNe where 
significant reverse shock interaction has likely not yet oc- 
curred, however, the simple model of purely toroidal magnetic 
field has been challenged by comparing X-ray observations of 
PWNe with the predictions of the Kennel & Coroniti (1984b) 
model. Chandra observations of 3C 58 (Slane et al. 2004) al- 
lowed comparison of the radial X-ray spectral index variation 
with the predictions of Kennel & Coroniti (1984b). For any 
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reasonable values of a these authors found the observed pho- 
ton index increased much slower in 3C 58 than expected from 
the model of Kennel & Coroniti (1984b). A similar deviation 
is observed for G21.5 - 0.9 (Slane et al. 2000), implying that 
the magnetic field structure deviates from a purely toroidal 
structure. 

The characteristic time scale for particle diffusion is de- 
termined by both the size of the nebula as well as diffusion 
coefficient D, such that r esc ps R 2 /6D. Alternately, the 
mean free path of a particle A = 3D/c. The electron dif- 
fusion coefficient is frequently parameterized in power-law 
form D{E e ) = D (E e /E efi ) s . For Bohm diffusion 5=1, 
while energy independent diffusion of course corresponds to 
5 = 0. To account for diffusion in a non-toroidal magnetic 
field we select a Bohm-type diffusion with 6=1, such that 
mean free path and escape time scale as: 

A = a E e 

T esc = R 2 {2acE e )- 1 (22) 

with a left as a free parameter in the model such that the mean 
free path is allowed to vary from the Bohm value to account 
for complex nebular field structure. For simplicity we ignore 
any perturbation the magnetic field might have on the escape 
timescale, since clearly particles must primarily move along 
rather than across magnetic field lines. 

3.8. Photon Fields 

Common practice in PWN SED modeling has been to ad- 
just the dust IR energy density to maximize the agreement 
between model and data. We instead fix the photon densi- 
ties at published values, reducing the number of model fit 
parameters. Adopting three primary photon fields (CMBR, 
far IR (dust), and starlight), we employ the interstellar radi- 
ation mapcube within the GALPROP suite (Porter & Strong 
2005) to estimate the photon fields at the Galactic radius of 
PSR J1826-1334. A distance of 4 kpc in the direction of 
PSR J 1826— 1334 corresponds to a Galactic radius of 4.7 kpc; 
at this radius, dust IR photons typically peak at T ps 32 K 
with a density of ss 1 eV cm~ 3 , while stellar photons peak at 
T ps 2500 K with a density of pa 4 eV cm -3 . At these den- 
sities IR and CMB photons contribute approximately equally 
to TeV inverse Compton emission, while stellar photons are 
inconsequential. Synchrotron self-Compton emission is neg- 
ligible for the magnetic fields present in HESS J1825— 137. 

3.9. SED Model 

A number of recent papers have applied various models to 
investigate the broadband emission from PWNe (Lemiere et 
al. (2009), Gelfand et al. (2009), Bucciantini et al. (2010), 
Fang & Zhang (2010), Slane et al. (2010), Tanaka & Taka- 
hara (2010)). Gelfand et al. (2009), in particular, constructed 
a detailed 1-D dynamical model of PWNe magnetic field, 
pressure, size, energetics and spectral evolution inside a host 
SNR. These authors demonstrated their model by fixing pa- 
rameters to mimic the Crab Nebula. While this model cer- 
tainly helps shed light on compact PWNe such as the Crab, the 
markedly different sizes of X-ray and TeV emission regions in 
HESS J1825-137 (along with most other middle aged Vela- 
like PWN) pose a challenge to 1-D models. 

We therefore treat the nebula as a series of twelve expand- 
ing bubbles in order to match the twelve H.E.S.S. spatial re- 
gions. Each of these bubbles is treated as homogeneous, and 




FIG. 5. — Snapshots of the evolution of nebular bubbles over three injection 
epochs. Dark gray corresponds to injection from the pulsar, light gray denotes 
the cooling phase, and arrows indicate diffusing particles. 



filled with particles (assumed to be electrons) as well as mag- 
netic field. We compute SEDs from evolving the electron 
populations over various lifetimes in a series of time steps, ex- 
tending the single-zone approach taken in (Abdo et al. 2010), 
(Ackermann et al. 201 1), (Grondin et al. 201 1). 

The twelve-zone modeling approach permits investigation 
of the twelve H.E.S.S. spectral extraction regions, yet the 
number of spatial zones need not be limited to the number 
of data regions and in theory one could utilize any number 
of zones and then project those zones onto the sky for SED 
construction. The spatial grid size is completely arbitrary, 
and hence the final SED plot should be independent of grid 
size. A finer grid of zones vastly increases computation time, 
however, and provides little new information given the limited 
number of data regions. Nevertheless, as a consistency check 
we double the number of zones (to 24) and plot the resultant 
SED. As expected, the model curves change imperceptibly in 
all but the innermost region. The slight alteration in Region 1 
is due to the steep variation in the magnetic field near the pul- 
sar. This lack of dependence of the SED on the spatial zone 
size provides a check of the SED model. 

3.10. Model Implementation 

Incorporating diffusion and projection effects between 
zones necessitates that all zones evolve in parallel, rather than 
independently. Evolution begins with Zone 12ati = t il2 = 
0, with escaping particles simply lost from the nebula. Once 
injection ends for Zone 12, at ti \i, Zone 1 1 begins to evolve, 
with its escaping particles injected into Zone 12. Once the 
injection period for Zone 11 ends, Zone 10 begins the in- 
jection phase, with Zones 11 and 12 in the cooling phases. 
This process is repeated until the current age of the pulsar is 
reached, with all twelve zones evolving simultaneously. Fig- 
ure 5 shows a schematic of the evolution of the nebula for 
three spatial zones. 

To compute the PWN SED, injected electrons are evolved 
through both the injection and cooling phases of each bub- 
ble's lifetime, with the age of each phase computed via Equa- 
tion 13. Each zone shares the same injection spectrum, mod- 
ulated by the spin down power. Cooling is computed from 
synchrotron (Blumenthal & Gould 1970, Equation 4.5), in- 
verse Compton (Blumenthal & Gould 1970, Equation 2.56), 
and adiabatic losses (Equation 17). Evolution occurs in a se- 
ries of time steps, such that during the injection phase at each 
time step we perform the following procedures: 

1. Calculate the size of the zone by integrating Equation 
16 

2. Compute the mean magnetic field by integrating over 
the volume of the zone via Equation 20 
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3. Inject pulsar wind particles into the zone according to 
Equation 10 

4. Inject particles diffusing from adjacent zones according 
to Equation 2 1 (if including diffusion) 

5. At each particle energy E, compute SE from syn- 
chrotron, inverse Compton, and adiabatic losses 

6. Calculate the subsequent particle spectrum at time t + 

St: E(t + St) = E - 5E 

7. Remove escaping particles according to Equation 21 (if 
including diffusion) 

Once injection ends, cooling is computed for the remainder 
of the evolution, following the same steps above (save for step 
3). At each time step this process is repeated for all zones with 
ti > t. When evolution completes, the synchrotron (Longair 
2010, Equation 8.58), and IC (Blumenthal & Gould 1970, 
Equation 2.48) fluxes are computed from the final electron 
spectrum and magnetic field for each zone. Finally, the flux 
within each region is tabulated using the coefficients of Table 
5. 

Model fitting is achieved by minimizing the \ 2 between 
model and data using the downhill simplex method described 
in Press et al. (1992). Data consists of both X-ray and 
H.E.S.S. spectral flux points (77 in total): 12 in Region 1, 
1 1 in Region 2, 9 in Region 3, and 5 data points for all exte- 
rior regions. For each ensemble of N variable parameters we 
evolve the system over the pulsar lifetime and calculate x 2 be- 
tween model curves and flux data points. The simplex routine 
subsequently varies the parameters of interest to minimize the 
fit statistic. We estimate parameter errors by computing \ 2 for 
a sampling of points near the best fit values and using these 
points to fit the N— dimensional ellipsoid describing the sur- 
face of A\ 2 = 2.71. Under the assumption of Gaussian er- 
rors, the projected size of this A\ 2 = 2.71 ellipsoid onto each 
parameter axis defines the 90% multi-parameter (projected) 
error. 

4. RESULTS 

4. 1 . Uniform Expansion, Constant Magnetic Field ( Model 1 ) 

We begin by considering a very simple model of con- 
stant magnetic field (~ 3 fiG), uniform expansion (v(r, t) = 
constant), and ignore diffusive effects. Uniform expansion 
implies that each zone undergoes expansion for 1 / 12 th of the 
pulsar lifetime. Injection begins at pulsar birth (t = 0) for 
Zone 12, while Zone 1 begins evolution at T — T/12, where 
T is the current age of the pulsar. The large size of the PWN 
implies a relatively high age, and the braking index and Pq se- 
lected in Table 6 give a total age of 21 kyr. This corresponds 
to a mean expansion velocity of 3900 km s -1 . We initially 
adopt a minimum particle energy E m i n = 100 GeV, and elec- 
tron injection fraction r\ = 0.8, which is a reasonable value 
for the expected particle dominated wind. 

The model gives a poor match to the data, apparent in the 
SED plot of Figure 6. Applying the full fitting machinery to 
five variables (magnetic field B, initial spin period Po, brak- 
ing index n, electron power-law index p, electron cutoff E cut ) 
yields an extremely low initial spin period < 1 ms. Figure 6 
therefore fixes the braking index (n = 3) and initial spin pe- 
riod (Po = 15 ms) to reasonable values; the other three vari- 
ables are fit and 90% multi-parameter error estimates are com- 
puted, as shown in Table 6. The fit is still quite poor, and the 
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FIG. 6. — Model 1. Fitted simple spectral energy distribution model of the 
nebula with a constant magnetic field (2.9 fiG) throughout the nebula for all 
time, constant expansion velocity such that each zone undergoes injection for 
l/12 th of the pulsar age (21 kyr), and x 2 /d.o.f. = 575/74. The SED plot 
is split into two panels for clarity, with a color scale evolving from dark red 
(Region 1, youngest), to green, to blue, to violet (Region 12, oldest). The 
line type evolves as well, from solid (inner, youngest) to long dot-dashed, 
short dot-dashed, long dashed, short dashed, and finally dotted (outer, oldest). 
The upper panel shows the model curves and H.E.S.S. data for the inner six 
regions, and Suzaku data for the inner three regions, while the lower panel 
shows the outer six regions. The fit is clearly quite poor, requiring fewer 
particles in the inner regions (and a higher magnetic field), with less cooling 
and more power in the outer regions. The gradual decrease of the cooling 
break energy with increasing region is primarily due to energy dependent 
synchrotron losses. 

fit value of E cut = 3200 TeV is implausibly large. Inspection 
of Figure 6 indicates that the normalization of each individual 
zone is highly skewed. The innermost zone has far too much 
VHE flux, while Zone 12 at the nebular frontier has far too 
little VHE flux. 

4.2. Velocity Profile, Constant Magnetic Field (Model 2) 

Though each region occupies the same 6' width, the nor- 
malization issues of Figure 6 intimate that the injection time 
of the zones must vary with distance from the pulsar, and 
hence the flow speed in the nebula must not be constant. We 
therefore adopt a velocity profile in the nebula of v(r,t) = 
r a C(t) (Equation 16). The best fit with this model and six 
variables is shown in Table 7. The radial velocity profile re- 
sults in a lengthening of the injection phase with increasing 
radial distance from the pulsar, and the values n and Po yield 
an age of 28 kyr. This model is shown in Figure 7, is far su- 
perior to Figure 6, and reproduces the VHE steepening trend 
moderately well for the inner six regions. Yet the outer model 
curves are all far steeper than the VHE data, implying an ex- 
cess of synchrotron cooling for the ~ 40 TeV particles re- 
quired to upscatter CMB photons to the ~ 5 TeV highest en- 
ergy H.E.S.S. data point (Equation 4). As a result, in the out- 
ermost regions (presumably corresponding to the oldest elec- 
trons), electrons of energy ~ 40 TeV have almost completely 
burned off. In the innermost region the X-ray flux is signifi- 
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FIG. 7. — Model 2. Fitted SED model with a constant 5.6 fj,G magnetic field 
throughout the nebula for all time, v(r) oc r~ °- 36 such that the length of the 
injection phase increases with zone number, and x 2 /d.o.f. = 236/71. The 
model curve line scheme is described in Figure 6. The radial velocity profile 
improves the fit over Figure 6, though the poor fit to the innermost X-ray 
data and excessive cooling in the outer zones suggest a radial magnetic field 
profile may provide a better fit. 

cantly under-produced, suggesting the field must be > 5 /iG 
near the pulsar. Evidently a non-uniform magnetic field pro- 
file is required, with higher field in the interior to match X-ray 
fluxes and lower field in the exterior in order to reduce syn- 
chrotron cooling for the oldest electrons. 

In order to better visualize the results of the SED fitting, 
we also show the ratios between the model zone fluxes and 
spectral indices and the measured data values (Figure 8). The 
model flux is computed by integrating the flux in the desired 
energy band (either 1 — 8 keV or 0.25 — 10 TeV) and model 
indices are computed by fitting to a power-law in the appro- 
priate band. The model fluxes and indices are compared to 
the data from Tables 3 and 4. 

4.3. Velocity Profile, Magnetic Field Profile (Model 3) 
In hopes of improving the fit of Figure 7, we invoke a ra- 
dial magnetic field profile of the form B(r, t) oc E^) 1 / 2 
(Equation 20), such that as particles traverse the nebula they 
experience an ever decreasing magnetic field. Similar to sec- 
tion 4.2, we select six fit parameters, though instead of fitting 
the final magnetic field value we fit the magnetic profile index 
f3. The best fit model yields an age of 34 kyr, and Table 6 re- 
veals an improved fit statistic over the constant magnetic field 
model. The non-constant magnetic field does a much bet- 
ter job of matching the innermost X-ray data (see Figures 9 
and 10 ), though we still have difficulty matching the VHE 
indices in the outer (purple) zones. 



FIG. 8. — Model 2. Top: Ratio of model flux to the data listed in Tables 3 
and 4, with X-ray 1 — 8 keV flux ratios for the inner three zones (dotted) and 
VHE 0.25 — 10 TeV flux ratios (solid) for all twelve zones. Bottom: Power- 
law indices over the same energy bands, with the model steadily worsening 
in the the outer zones. 
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Energy (eV) 

FIG. 9.— Model 3. Fitted SED model with B(r, t) oc r" ' 74 £(r.) 1/2 , 
v(r) oc r~ 0- 60 , and \ 2 /d.o.f. = 217/71. The model curve line scheme is 
described in Figure 6. The radial magnetic field profile results in a superior 
fit to the innermost X-ray data compared to Figure 7, though we still suffer 
from excessive cooling losses in the outermost zones. 



12 



Van Etten & Romani 




Zone 



FIG. 10. — Model 3. Top: Ratio of model flux to the data listed in Tables 3 
and 4, with X-ray 1 — 8 keV flux ratios for the inner three zones (dotted) and 
VHE 0.25 — 10 TeV flux ratios (solid) for all twelve zones. The X-ray flux 
ratios are much improved over Figure 8. Bottom: Power-law indices, with 
the model steadily worsening in the the outer zones. 

4.4. Velocity Profile, Magnetic Field Profile, Diffusion 
(Model 4) 

The severe under-prediction of VHE flux and over- 
prediction of VHE index in the outer two zones of Figures 9 
and 10 indicates a general lack of electrons with E > 10 
TeV. One possibility of addressing this issue is simply reduced 
cooling by lowering the magnetic field. Yet even in Figures 7 
and 8 with a low 6 /iG field the Region 12 model is far steeper 
than the data due to excessive cooling. In addition, Figures 9 
and 10 demonstrate that higher fields are needed in the in- 
ner zones, which only serves to hasten synchrotron burn-off. 
Particle re-acceleration as particles traverse the nebula might 
explain the VHE flux ~ 80 pc from the pulsar, perhaps due to 
reverse shock interactions. Such acceleration is highly specu- 
lative, however, and beyond the scope of this paper. A more 
standard method of maintaining the high energy particle pop- 
ulation far from the pulsar, which we explore below, is rapid 
energy-dependent diffusion with mean free path A oc E e . 

For this model we select seven variables: velocity radial 
profile index a, magnetic profile index f3, initial spin period 
Pq, braking index n, electron power-law index p, electron ex- 
ponential cutoff E cut , and mean free path A. Figure 1 1 shows 
the SED of the best fit model, scaled to display the aggregate 
gamma-ray data. This aggregate data is discussed in Section 
5. Figure 12 displays the flux and index ratios for this model. 
The overall fit to the twelve spectral regions is far superior to 
previous models, particularly in the outer nebula. The largest 
departure of the model from the data is the X-ray spectral in- 
dex of Region 3, though we remind the reader that the very 
hard Region 3 X-ray index reported in Table 4 may reflect 
residual background contamination. Table 6 lists the fit pa- 
rameters, while Table 7 lists derived parameters from the fit 
(injection time, final magnetic field, and x 2 ). As expected 
from the scaling relations of section 3.1, the particle residence 



time of Zone 1 is well under 1000 years. The braking index 
n and initial spin period Pq yield an age of 40 ± 9 kyr and a 
mean expansion velocity of 2000 km s _1 . The mean free path 
of A = 1.8(i? e /100 TeV) parsec predicts an escape timescale 
ofr esc w 90(i?/10pc) 2 (E e /100TeV)- 1 year. This model 
does have some difficulty matching the H.E.S.S. flux and 
slope in Region 1, which is due to the rapid diffusive loss 
of high energy particles. Most likely, the simple diffusive 
model employed vastly oversimplifies the complex nature of 
this source, and diffusion coefficients vary in both space and 
time. Nevertheless, the x 2 /d.o.f. for Model 4 indicates that 
the general fit is quite good. 

As a check of our magnetic field model, we also allow 
the normalization at the termination shock to vary from its 
adopted value of 400 fiG. The best fit model with this extra fit 
parameter yields a termination shock field of 320 fiG with no 
improvement of the fit statistic (x 2 / d.o.f. = 118/69, versus 
X 2 /d.o.f . = 119/70 for the field fixed. Fitting the magnetic 
field normalization therefore does not yield a markedly better 
fit, and the termination shock value of 400 fiG appears quite 
reasonable. One further check of the magnetic field model 
is achieved by fitting a model with diffusion over uniform 
magnetic field (simply replacing the magnetic field radial in- 
dex of Model 4 with a constant magnetic field value). This 
model cannot adequately match the data, yielding a best fit 
of x 2 /d.o.f. — 429/70, and rules out uniform magnetic field 
as a viable option. The adopted IR photon density of ps 1 
eV cm' 3 could easily vary by a factor of ~ 2, though varying 
the photon density does not qualitatively change the results 
presented here. A higher (lower) IR photon density means 
fewer (more) high energy electron are required to match the 
HESS data, implying a slightly softer (harder) electron spec- 
tral index or lower (higher) high energy cutoff. 

4.4. 1 . Relativistic Maxwellian Plus Power-Law Injection Spectrum 

The SED model described in the previous section takes into 
account only X-ray and VHE H.E.S.S. data, though ideally 
the model should be consistent with Fermi LAT data as well. 
Note that while we do not fit to the LAT points in Figure 1 1 , 
the rough agreement between these data and the integrated 
(total) model flux in the LAT band is encouraging. However, 
the Fermi LAT photon index of T = 1.4 ± 0.2 (Grondin et al. 
201 1) predicts an electron index of p = 2T - 1 = 1.8 ± 0.4. 
This is somewhat harder than the electron index of p = 2.2 
utilized in Figures 1 1 and 12, and an injection spectrum with 
fewer low energy electrons might plausibly improve the fit to 
the LAT data. 

One option is adopting the relativistic Maxwellian plus 
power-law tail electron spectrum proposed by Spitkovsky 
(2008), and described in modeling terms by Grondin et al. 
(2011). Fitting with this model requires one additional pa- 
rameter (the temperature of the Maxwellian), bringing the to- 
tal up to eight. The best fit model with this injection spec- 
trum yields a fit somewhat worse than the power-law case 
(X 2 /d.o.f. = 123/69, versus 122/70 for the power-law spec- 
trum). In addition, an extremely short initial spin period of 
Pq ~ 2 ms is required, and the low kT = 0.2 TeV results 
in a severe over-prediction of GeV flux and a far worse fit to 
the LAT data than the power-law case. Grondin et al. (201 1) 
found an adequate fit to the aggregate data with a one-zone 
model, yet the inclusion of multiple zones seems to wash out 
any advantage of a relativistic Maxwellian electron spectrum. 
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TABLE 6 

MULTI WAVELENGTH SED FIT TO ALL REGIONS 



Model 


a (v oc r a ) 


B/P (B oc r?) 


Po (ms) 


n (Q oc Q n ) 


V 


Scut (TeV) 


^ ( 100 TcV ) pC 


X 2 /d.o.f. 


1 




2.9±0.3 a 


15 b 


3 b 


2.42 ± 0.06 


3800 ± 800 




575/74 


2 


-0.36 ±0.03 


5.6±0.4 a 


12 ±3 


2.5 ±0.1 


2.30 ±0.01 


130 ± 10 




236/71 


3 


-0.60 ±0.08 


-0.74 ±0.02 


14 ±7 


2.2 ±0.1 


2.31 ±0.02 


190 ± 30 




217/71 


4 


-0.51 ±0.09 


-0.69 ± 0.01 


13 ±7 


1.9 ±0.2 


2.24 ±0.02 


230 ± 90 


1.8 ±0.3 


119/70 



a Final (uniform) magnetic field (fiG) 
b Held fixed 




FIG. 1 1 . — Model 4. Fitted SED model, scaled to show the aggregate data, 
with B(r,t) Kr" - 69 ^) 1 / 2 .^) oc r" - 51 , A = 1.8 (E e /100 TeV) 
parsec, andx 2 /rf-o./. = 119/70. The model curve line scheme is described 
in Figure 6. The inclusion of rapid diffusion among zones provides enough 
high energy electrons in the outer nebula to match the H.E.S.S. data in Zones 
10 — 12. The slate data points denote the aggregate Fermi LAT and H.E.S.S. 
data, with both statistical and systematic error estimates. The slate dotted 
line (representing the sum of flux from all twelve zones) is scaled up by 25% 
to account for the greater total flux in the aggregate extraction region (see 
Section 2.2). The overall fit to the aggregate data with this scaling is quite 
impressive. 



TABLE 7 
Derived Parameters for Model 4 



Zone 


R (pc) 


Sti (year) a 


B ( M G)» 


x 2 


# data points 


1 


7 


600 


12 


23 


12 


2 


11 


1200 


6.6 


13 


11 


3 


21 


1500 


4.8 


6.6 


9 


1 


28 


1900 


3.8 


1.6 


5 


5 


35 


2200 


3.2 


8.0 


5 


6 


42 


2500 


2.8 


10 


5 


7 


49 


2900 


2.5 


8.9 


5 


8 


56 


3300 


2.2 


11 


5 


9 


63 


3800 


2.1 


7.4 


5 


10 


70 


4400 


1.9 


6.1 


5 


11 


77 


5500 


1.8 


16 


5 


12 


84 


9900 


1.7 


4.4 


5 



a Injection timescale 
b Final magnetic field 



FIG. 12. — Model 4. Top: Ratio of model flux to the data listed in Tables 
3 and 4, with X-ray 1 — 8 keV flux ratios for the inner three zones (dotted) 
and VHE 0.25 — 10 TeV flux ratios (solid) for all twelve zones. Bottom: 
Power-law indices. For both upper and lower panels the outer-zone ratios are 
far superior to Figures 8 and 10. 

5. DISCUSSION 

As an initial check of the validity of the scalings applied to 
both the Suzaku and H.E.S.S. data, we refit Model 4 with the 
raw data. Using the original H.E.S.S. data points and X-ray 
data taken from regions coinciding with the HESS quarter an- 
nuli, the best fit parameters are within lcr of the values listed 
in Table 6 for Model 4, except for the magnetic field index. 
The different H.E.S.S. flux ratios between interior zones re- 
sults in a magnetic field index that differs from the Model 4 
value by +0.02, or 2a. Despite the similarity in fit param- 
eters, however, the fit statistic increases dramatically, from 
X 2 /d.o.f. = 119/70 for Model 4, to x 2 Mo./. = 169/70. 
Nearly all of this increase in x 2 comes from the inner two re- 
gions. The H.E.S.S. excess map (Figure lb) reveals that there 
appears to be significantly more flux 6' — 12' from the pul- 
sar than within 6 of the pulsar. Yet due to the small original 
H.E.S.S. extraction region for the 6' — 12' region (which is 
25% smaller than the — 6' region), the reported H.E.S.S. 
fluxes are identical for the inner two regions. The use of iden- 
tical fluxes for the inner two zones therefore results in a worse 
fit for the model, which assumes radial variations in the flow 
speed and a greater number of particles farther from the pul- 
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sar. Nevertheless, the commonality of the fit parameters with 
the raw and scaled data provides some support for the scaled 
approach taken throughout this paper. One further check of 
the X-ray scaling is to fit the scaled H.E.S.S. data and raw, 
unsealed X-ray data. With this slightly lower X-ray flux the 
best fit parameters are well within errors for Model 4, and give 
a nearly identical % 2 . 

Of the four models scenarios described above, the final 
model which incorporates diffusion provides by the far the 
best fit to the data. The relativistic Maxwellian injection spec- 
trum yields a somewhat worse fit statistic to the simple power- 
law case, and though the power-law injection spectrum also 
better matches the Fermi LAT data, definitive discrimination 
between lepton injection spectra is not possible at this stage. 
Yet we can convincingly argue that the pulsar wind must in- 
ject significant numbers of electrons up to energies of at least 
100 TeV, the low energy power-law index must be p ~ 2 in 
order to match the Fermi LAT data, that a high percentage 
of the pulsar spin-down energy must go into leptons, and that 
some mechanism (we invoke diffusion) maintains high energy 
particles far from the pulsar. In the following discussion we 
consider the simple power-law spectrum injection spectrum 
(Model 4) shown in Figures 1 1 and 12. 

The sum of fluxes from the twelve spatial zones has been 
scaled up some 25% in Figure 1 1 due to the fact that our spec- 
tral extraction regions do not capture the entirety of the counts 
within the aggregate nebula. This scaled flux ideally should 
match the Fermi LAT and H.E.S.S. aggregate data points 
(which both represent larger spatial regions), which appears 
to be the case for Figure 1 1 . Computing the fit statistic for the 
aggregate H.E.S.S. and LAT data gives x 2 = 27 for the 16 
gamma-ray data points (assuming only statistical errors on the 
data points, and ignoring the LAT upper limit). Given the very 
small statistical errors on the H.E.S.S. data points the agree- 
ment between gamma-ray data and the summed model curves 
is remarkable. For comparison, the best fit model utilizing the 
relativistic Maxwellian injection spectrum nets a \ 2 — 50 for 
the 16 gamma-ray data points. The 'missing' « 25% of flux 
implies a missing population of particles contributing to the 
IC flux of the nebula. Such particles have likely diffused to 
the north and east of the pulsar exterior to our extraction re- 
gions. A tweak of the initial pulsar spin-down from 13 ms to 
1 1 ms (well within errors) provides w 40% greater spin-down 
energy, more than enough to enough to raise the model bulk 
nebular flux to the observed level. 

The models proposed above utilize r\ = 0.8, and we find an 
adequate fit to the X-ray and H.E.S.S. data for P — 13 ms. 
Neglecting the 'missing' energy component, this connotes a 
total of 1.1 x 10 50 erg of pulsar spin-down energy. Disregard- 
ing losses, we therefore expect of order 10 49 erg in magnetic 
field energy. While the magnetic field evolution and profile 
adopted in our model does not take into account the total en- 
ergy budget of the pulsar, a reasonable model should never- 
theless contain <~ 10 49 erg of magnetic energy for the choice 
of r\ = 0.8. Computing the magnetic field energy from the fi- 
nal magnetic field profile B(r, T) = 400(r/0.03 pc)" 69 ^G 
yields an energy of 4.4 x 10 48 erg, which is quite reason- 
able given that the magnetic field likely experiences some 
losses due to expansion. Even though the vast majority of the 
pulsar energy initially goes into particles, summing the en- 
ergy content of the final electron spectrum reveals that elec- 
tron cooling losses leave only 6.9 x 10 48 erg remaining in 
the form of leptons. We can also compute the value of a 



from the assumed termination shock radius of 0.03 pc and 
400 [iG field strength. Adopting a flow speed at the termi- 
nation shock of c/3 and a spherical geometry, one achieves 
an estimate of 7 x 10 35 ergs _1 of magnetic energy flux at 
the termination shock. For the pulsar spindown value of 
E = 2.8 x 10 36 ergs -1 this yields a w 0.3. This value is 
consistent with the adopted value of r] = 0.8, which corre- 
sponds to a value of a = (1 — r\)jr\ = 0.25. One further 
consistency check for our model comprises comparing the 
available spin-down power to the total energy present in the 
nebula plus cooling and escape losses. Indeed, summing the 
leptonic energy, magnetic energy, cooling losses, and escape 
losses yields 1.1 x 10 50 erg, precisely the energy available 
from the pulsar. 

The energy lost due to adiabatic cooling goes into acceler- 
ating a shell of swept-up material, and we can use estimates of 
this swept-up mass as another check of the validity of the SED 
model. For the adopted SNR density profile of a constant den- 
sity inner core surrounded by a p oc r~ 9 outer envelope the 
transition between the two density regimes propagates out at 
a constant velocity v t , such that (Blondin et al. 2001): 

Vt „ 6 500 ( - ( ^y 1/2 kms -i.(23) 

^3x10 51 erg/ \SM Q J 

This is much faster than the Model 4 outer PWN expansion of 
2000 km s -1 , so we can treat the PWN as always expanding 
into the constant density ejecta core, with density Blondin et 
al. (2001): 

x (jofe) 3 (24 > 

The amount of mass swept up by the PWN over its lifetime 
can be approximated as: 

M sw = J p ej (t)(4Trr 2 dr, (25) 

with ti the initial timestep, r — v outer t and ( the filling fac- 
tor of the PWN, equal to 0.31 for the adopted morphology. 
Integrating Equation 25 over the 40 kyr lifetime of the pul- 
sar provides a very rough estimate of 1M Q of swept up mass. 
Accelerating 1 Mq to 2000 km s -1 (the mean expansion ve- 
locity for a pulsar age of 40 kyr) requires some 4 x 10 49 erg. 
For the parameters selected in Model 4, the SED model posits 
a total adiabatic loss of 3.4 x 10 49 erg, quite close to the value 
quoted above. Recall that in the SED model adiabatic losses 
depend only on the assumed PWN size evolution (Equation 
17), and do not take into account any SNR properties. There- 
fore the fact that the two adiabatic loss estimates are not glar- 
ingly different provides some support to the soundness of the 
SED model. 

The computed pressure within the entire nebula can be 
used as a further sanity check on our model by comparing 
to the expected PWN pressure after 40 kyr. Given the de- 
gree of particle cooling, the model's final particle pressure 

(Ppwn = 3 y" = L0 x 10 -13 ergcm~ 3 ) is less than the 

magnetic pressure (J^ = 1.9 x 10~ 13 erg cm -3 ). Within a 
freely expanding SNR the PWN pressure scales as (van der 
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Swaluw et al. 2001): 

Ppwn ~ 3/25 p ej (t)v pwn , (26) 

with p e j (t) given by Equation 24, and v pwn the PWN expan- 
sion velocity. After 40 kyr this predicts a mean nebular pres- 
sure of Ppwn ~2x 10~ 13 erg cm~ 3 , in good agreement with 
the computed values. 

The adopted termination shock radius of 0.03 pc predicts 
a pressure of m 1 x 10 -9 erg cm -3 at the termination shock 
(Equation 19), significantly higher than the mean pressure of 
~ 3x 10 -13 erg cm -3 . This high pressure at the shockrapidly 
dissipates away from the pulsar, however, with the bulk of the 
nebula nearly isobaric at late times; each of the 10 outer zones 
has a pressure within ps 50% of the mean nebular pressure, 
as evidenced in the pressure and energy plot of Figure 13. 
This plot illustrates that the magnetic pressure exceeds parti- 
cle pressure for the inner regions, and that the total particle en- 
ergy exceeds the bulk magnetic field energy in the nebula. The 
total magnetic field pressure scales as Ppwn b oc t 3 , while 
the total particle pressure evolves as P pwn , P oc t~ 2 9 . This is 
precisely the P pwn oc t~ 3 evolution predicted by Equations 
24 and 26 for a freely expanding PWN. 

In addition to the PWN pressure and energy we can also 
compute a variety of quantities which vary over the pulsar 
lifetime. We calculate the mean magnetic field evolution as 
B oc t^ 1 - 6 , which is quite similar to the B oc evolution 
computed by Gelfand et al. (2009) in their single-zone model 
during the free expansion phase. The various luminosities in 
the nebula are plotted in Figure 14, revealing that at late times 
the total luminosity of the nebula exceeds the pulsar spin- 
down power. This figure also demonstrates that after a few 
thousand years adiabatic losses overtake synchrotron radia- 
tion to provide the dominant cooling mechanism, though par- 
ticles of energy E e > 10 TeV primarily lose energy from syn- 
chrotron and IC cooling since for both the cooling timescale 
r oc E~ 1 . Indeed, at late time the IC cooling actually exceeds 
the synchrotron cooling due to the low mean magnetic field of 

We construct radial flux profiles by summing the model 
fluxes in the Suzaku and H.E.S.S. wavebands within each re- 
gion, which can be compared with the data points plotted in 
Figure 2. We are also free to compute the profile in any arbi- 
trary energy range, and so in Figure 15 we plot the projected 
surface brightness in the Fermi LAT 1-100 GeV range. Since 
the LAT energy range lies in the uncooled regime, the flux pri- 
marily tracks the energy content of electrons. As a result, the 
LAT flux peaks toward the outer realm of the nebula where 
the majority of electrons were injected when the pulsar was 
much younger with a far higher E, and this gives that a rather 
uniform counts map in the LAT energy band. 

Comparison of point source and extended Gaussian mor- 
phologies in the 10 — 100 GeV range by Grondin et al. (201 1) 
yielded a significance of m 8.5a for a source extension of 
0.56° ± 0.07° centered 0.5° southwest of the pulsar. Similar 
results are achieved with a uniform disk morphology, though 
the extension increases to 0.67° ± 0.02°. Fitting the source 
to a spatial template provided by the H.E.S.S. excess map 
(E > 200 GeV) decreases the test statistic by 42 sa 6.5a, 
implying that the LAT emission shows a distinctly different 
morphology from the H.E.S.S. emission. Indeed, the LAT 
source extension of ~ 0.6° centered 0.5° southwest of the 
pulsar is remarkably consistent with Figure 15. 
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FIG. 13. — Top: Magnetic field pressure (dotted) and particle pressure 
(solid) of each zone. Total pressures are marked by the thick black lines. 
Bottom: Magnetic field energy (dotted) and particle energy (solid) of each 
zone. Black lines mark total energies. In both panels the values correspond 
to Model 4. 
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FIG. 14. — Bulk luminosity of the nebula for Model 4. The solid gold line 
denotes pulsar spin-down, while the solid black line shows the total lumi- 
nosity from all processes. Synchrotron (green long-dashed), IC (blue short- 
dashed), and adiabatic (red dotted) cooling is also shown. The violet line 
denotes diffusive escape losses from the nebula. 

6. CONCLUSIONS 

Analysis of recently public Suzaku X-ray data covering the 
majority of the very bright VHE source HESS J1825— 137 in- 
dicate that X-rays extend only ps 1/4 as far as VHE gamma- 
rays. The variable extent of the nebula in different wave- 
bands, as well as the observed VHE spectral softening away 
from the pulsar, is naturally explained by electron cooling 
losses. Time-dependent 3-D modeling of the twelve data re- 
gions allows one to constrain the electron injection proper- 
ties, spin-down behavior of the pulsar, nebular velocity pro- 
file, and magnetic field profile. A choice of 80% of the 
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FIG. 15. — Nebula surface brightness profile indicating a remarkably uni- 
form Fermi LAT profile. Red data points denote Suzaku X-ray data from 
Table 4 (statistical plus sytematic errors shown), while the red curve shows 
the Model 4 flux prediction in the 1-8 keV band. Blue data points mark the 
H.E.S.S. flux from Table 3, with the blue curve the 0.25-10 TeV model pre- 
diction. The model 1-100 GeV band surface brightness is shown in green. 

pulsar spin-down going into electrons is consistent with the 
data, though due to particle cooling at the current epoch mag- 
netic field energy is w 1/2 x the particle energy in the neb- 
ula. This cooling also results in a higher magnetic than par- 
ticle pressure, with a nearly isobaric particle pressure pro- 
file in the outer ten zones and the total pressure falling as 
Ppwn <x t~ 3 - For the best fit model (Model 4) an initial 
spin period of m 13 ms provides adequate pulsar spin-down 
energy to power the nebula, and combined with the braking 
index of n = 1.9 this predicts an age of 40 kyr and an expan- 
sion velocity of 2000 km s _1 . The interior velocity profile is 
found to scale as v(r) ~ r~ - 5 , while the magnetic field pro- 
file scales as B(r,t) oc r~ -° J E(tf/ 2 , with the field falling 



from m 400 /iG at the termination shock to m 2 [iG in the far 
reaches of the nebula and the mean field scaling as B oc t~ 16 . 

The Fermi LAT photon index of 1 .4 connotes an electron in- 
dex near the canonical 2.0. However, fitting to the X-ray/TeV 
and the higher energy portion of the electron spectrum, we 
find a somewhat softer power law index of p w 2.2. A 
relativistic Maxwellian electron spectrum provides a some- 
what worse fit to the X-ray and H.E.S.S. data compared to 
a simple power-law injection, and the power-law case re- 
quires fewer fit parameters and better matches the Fermi LAT 
data. Regardless of injection spectrum, the existence of 10 
TeV photons in the outer reaches of the nebula requires sub- 
stantial particle diffusion; adopting a diffusion timescale of 
r esc w 90(i?/10pc) 2 (E c /100TeV)- 1 year allows an ade- 
quate fit to the data. This is in contrast to the common assump- 
tion of a purely toroidal PWN magnetic field which could 
contain particles extremely efficiently. The model predicts a 
rather uniform Fermi LAT surface brightness out to ~ 1 from 
the pulsar, in good agreement with the recently discovered 
LAT source centered 0.5° southwest of PSR J1826— 1334, 
with extension 0.6 ± 0.1°. 

The multi-zone time-dependent SED model fitting ap- 
proach applied in the paper extends simple one-zone models 
to allow investigation of spatially dependent properties such 
as flow speed, magnetic field, photon index, and flux levels. 
The growing number of sources with spatially resolved X-ray 
and VHE measurements (e.g. Vela-X, HESS J1303-631) are 
therefore prime targets for such multi-zone modeling. 
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